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ABSTRACT 

We propose a deconvolution algorithm for images blurred 
and degraded by a Poisson noise. The algorithm uses a fast 
proximal backward-forward splitting iteration. This iter- 
ation minimizes an energy which combines a non-linear 
data fidelity term, adapted to Poisson noise, and a non- 
smooth sparsity-promoting regularization (e.g £i-norm) 
over the image representation coefficients in some dictio- 
nary of transforms (e.g. wavelets, curvelets). Our results 
on simulated microscopy images of neurons and cells are 
confronted to some state-of-the-art algorithms. They show 
that our approach is very competitive, and as expected, 
the importance of the non-linearity due to Poisson noise 
is more salient at low and medium intensities. Finally an 
experiment on real fluorescent confocal microscopy data 
is reported. 

Index Terms — Deconvolution, Poisson noise, Confo- 
cal microscopy, Iterative thresholding, Sparse representa- 
tions. 

1. INTRODUCTION 

Fluorescent microscopy suffers from two main sources of 
degradation: the optical system and the acquisition noise. 
The optical system has a finite resolution introducing a 
blur in the observation. This degradation is modeled as 
convolution with the point spread function (PSF). The sec- 
ond source of image degradation is due to Poisson count 
process (shot noise). In presence of Poisson noise, sev- 
eral deconvolution algorithms have been proposed such 
as the well-known Richardson-Lucy (RL) algorithm or 
Tikhonov-Miller inverse filter, to name a few. RL is ex- 
tensively used for its good adaptation to Poisson noise, 
but it tends to amplify noise after a few iterations. Reg- 
ularization can be introduced in order to avoid this issue. 
In biological imaging deconvolution, many kinds of reg- 
ularization have been suggested: total variation with RL 
[ 1 ] which gives staircase artifacts, Tikhonov with RL (see 
|2] for a review), etc. Wavelets have also been used as a 
regularization scheme when deconvolving biomedical im- 
ages; | 3 ] presents a version of RL combined with wavelets 
denoising, and |4] use the thresholded Landweber itera- 
tion of 1 5 1. The latter approach implicitly assumes that the 
contaminating noise is Gaussian. 

In the context of deconvolution with Gaussian white 
noise, sparsity-promoting regularization over orthogonal 



wavelet coefficients has been recently proposed |5|. 
Generalization to frames was proposed in (71 [8). In |9|, 
the authors presented an image deconvolution algorithm 
by iterative thresholding in an overcomplete dictionary 
of transforms. However, all sparsity-based approaches 
published so far have mainly focused on Gaussian noise. 

In this paper, we propose an image deconvolution al- 
gorithm for data blurred and contaminated by Poisson 
noise. The Poisson noise is handled properly by using the 
Anscombe variance stabilizing transform (VST), leading 
to a non-linear degradation equation with additive Gaus- 
sian noise, see (|2]). The deconvolution problem is then 
formulated as the minimization of a convex functional 
with a non-linear data-fidelity term reflecting the noise 
properties, and a non- smooth sparsity-promoting penalty 
over representation coefficients of the image to restore, 
e.g. wavelet or curvelet coefficients. Inspired by the work 
in |6|, a fast proximal iterative algorithm is proposed to 
solve the minimization problem. Experimental results are 
carried out to compare our approach on a set of simu- 
lated and real confocal microscopy images, and show the 
striking benefits gained from taking into account the Pois- 
son nature of the noise and the morphological structures 
involved in the image. 

Notation 

Let TL a real Hilbert space, here a finite dimensional vector subspace of 
W 1 . We denote by || . || 2 the norm associated with the inner product in 7i, 
and I is the identity operator on H. x and a are respectively reordered 
vectors of image samples and transform coefficients. A function / is co- 
ercive, if limna-u _> +00 / (x) = +oo. To(H) is the class of all proper 
lower semi-continuous convex functions from 7i to ] — oo, +oo]. 

2. PROBLEM STATEMENT 

Consider the image formation model where an input image 
x is blurred by a point spread function (PSF) h and con- 
taminated by Poisson noise. The observed image is then a 
discrete collection of counts y = (yi)i^i^ n where n is the 
number of pixels. Each count jji is a realization of an in- 
dependent Poisson random variable with a mean (h ® x)i, 
where © is the circular convolution operator. Formally, 
this writes yi ~ V {{h® x)i). 

A naive solution to this deconvolution problem would 
be to apply traditional approaches designed for Gaussian 
noise. But this would be awkward as (i) the noise tends 
to Gaussian only for large mean (h ® x)i (central limit 



theorem), and (ii) the noise variance depends on the mean 
anyway. A more adapted way would be to adopt a bayesian 
framework with an appropriate anti-log-likelihood score 
reflecting the Poisson statistics of the noise. Unfortunately, 
doing so, we would end-up with a functional which does 
not satisfy some key properties (the Lipschitzian prop- 
erty stated after ©), hence preventing us from using the 
backward-forward splitting proximal algorithm to solve 
the optimization problem. To circumvent this difficulty, 
we propose to handle the noise statistical properties by 
using the Anscombe VST defined as 

Zi = + 1 < i < n. (1) 

Some previous authors ifTOl have already suggested to use 
the Anscombe VST, and then deconvolve with wavelet- 
domain regularization as if the stabilized observation Z{ 
were linearly degraded by h and contaminated by additive 
Gaussian noise. But this is not valid as standard asymp- 
totic results of the Anscombe VST state that 

* i = 2^(fc®a;) i + § + e, e ~ AT(0,1) (2) 

where e is an additive white Gaussian noise of unit vari- 
anceQ. In words, z is non-linearly related to x. In Sec- 
tion |U we provide an elegant optimization problem and 
a fixed point algorithm taking into account such a non- 
linearity. 

3. SPARSE IMAGE REPRESENTATION 

Let x G H be an y/n x y/n image, x can be written as 
the superposition of elementary atoms <p 7 parametrized 
by 7 G X such that x = ^ ie x a i i Pi ~ l-^l = 

L, L ^ n. We denote by $ the dictionary i.e. the n x 
L matrix whose columns are the generating waveforms 
((/? 7 ) 7GX all normalized to a unit ^2-norm. The forward 
transform is then defined by a non-necessarily square ma- 
trix T = $ T G R Lxn . In the rest of the paper, $ will be 
an orthobasis or a tight frame with constant A. 

4. SPARSE ITERATIVE DECONVOLUTION 

4.1. Optimization problem 

The class of minimization problems we are interested in 
can be stated in the general form [6) : 

arg min fx (x) + f 2 (x) . (3) 

xen 

where /i G T (H), h € r (W) and /i is differentiable 
with ft-Lipschitz gradient. We denote by M the set of 
solutions. 

From 0, we immediately deduce the data fidelity term 

FoHo$ (a), with (4) 



Vigorously speaking, the equation is to be understood in an asymp- 
totic sense. 



where H denotes the convolution operator. From a statisti- 
cal perspective, dU) corresponds to the anti-log-likelihood 
score. 

Adopting a bayesian framework and using a standard 
maximum a posteriori (MAP) rule, our goal is to minimize 
the following functional with respect to the representation 
coefficients a\ 

(Pa,vO : ar & min J(a) (5) 

a 

L 

J : a h-> F o H o $ (a) + %c ° $ (a) + A i/>(<Xi), 

V v ' 

/2(a) 

where we implicitly assumed that (^)i^^l are indepen- 
dent and identically distributed. The penalty function ip 
is chosen to enforce sparsity, A > is a regularization 
parameter and %c is the indicator function of a convex set 
C. In our case, C is the positive orthant. We remind that 
the positivity constraint is because we are fitting Poisson 
intensities, which are positive by nature. 

4.2. Proximal iteration 

We now present our main proximal iterative algorithm to 
solve the minimization problem (Pa,v>) : 

Theorem 1. (Pa,^) nas at least one solution (M ^ 0). 
The solution is unique if ip is strictly convex or if & is a 
orthobasis and Kei(H) = 0. For t ^ 0, let (nt)t be sucn 
thatO < inf t < sup t pi t < (§) 3/2 I (2A \\H\\ 2 2 \\z\Q. 
Fix ao G C o <l>, for every t ^ 0, set 

ott+i = prox Mt/2 (a t - inVhipLt)) , (6) 

where V/i is the gradient of f\ and prox^^ is computed 
using the following iteration: let v t (l — v t ) = +oo, 
take 7° G W, and define the sequence of iterates: 

7 £+1 = f + vt (rprax^^i || _ a||2 orprax. c/ -i) ( 7 *), 

w/^prox^^ _ a||2 ( 7 *) = ^prox^A^((^+7|)/2) 

= prox 2c/ = A _1 ^> T o? c oH(l - A' 1 ^ o $), 
rprox^ = 2 prox^ —I <md is the projector onto the 
positive orthant (Vcv)i ~ max (7^, 0). Then, 

7 £ 7 and prox /it/2 (a) = TV (7). (8) 

Then (a t )t^o converges (weakly) to a solution of(P\^). 

A proof can be found in [ 1 1 ]. prox 5 ^ is given by, 

Theorem 2. Suppose that ( i) ip is convex even- symmetric , 
non-negative and non-decreasing on [0, +00), and ip(0) — 
0. (ii) ip is twice differentiable on R \ {0}. (Hi) ip is con- 
tinuous on R, it is not necessarily smooth at zero and ad- 
mits a positive right derivative at zero ip+(0) > 0. Then, 



the proximity operator prox^^) = a(/3) has exactly one 
continuous solution decoupled in each coordinate : 

a(6) = l° mi<^(0) (9) 

lK J [Pi-Srlti&i) if >^+(0) 

See 0. Among the most popular penalty functions 
i/j satisfying the above requirements, we have ip(ai) = 
\cti\, in which case the associated proximity operator is 
soft- thresholding. Therefore, © is essentially an iterative 
thresholding algorithm with a positivity constraint. 

5. RESULTS 

The performance of our approach has been assessed on 
several datasets of biological images: a neuron phantom 
and a cell. Our algorithm was compared to RL with 
total variation regularization (RL-TV [1]), RL with multi- 
resolution support wavelet regularization (RL-MRS L12J), 
the naive proximal method that would assume the noise 
to be Gaussian (NaiveGauss L4j), and the approach of 
|[T0l (AnsGauss). For all results presented, each algorithm 
was run with 200 iterations, enough to reach convergence. 
Simulations were carried out with an approximated but re- 
alistic PSF [ 13 ] whose parameters are obtained from a real 
confocal microscope. As usual, the choice of A is crucial 
to balance between regularization and deconvolution. For 
all the situations below, A was adjusted in order to reach 
the best compromise. 

In FigHJa), a phantom of a neuron with a mushroom- 
shaped spine is depicted. The maximum intensity is 30. Its 
blurred and blurred+noisy versions are in (b) and (c). With 
this neuron, and for NaiveGauss, AnsGauss and our ap- 
proach, the dictionary contained the wavelet orthogonal 
basis. The deconvolution results are shown in FigQId)- 
(h). As expected the worst results are for the AnsGauss 
version, as it does not take care of the non-linearity of the 
Anscombe VST. RL-TV shows rather good results but the 
background is full of artifacts. Our approach provides a vi- 
sually pleasant deconvolution result on this example. It ef- 
ficiently restores the spine, although the background is not 
fully cleaned. RL-MRS also exhibits good deconvolution 
results. Qualitative visual results are confirmed by quan- 
titative measures of the quality of deconvolution, where 
we used both the i\ -error (adapted to Poisson noise), and 
the traditional MSE criteria. The l\ -errors for this image 
are shown by Tab. Q] (similar results were obtained for the 
MSE). In general, our approach performs well. It is com- 
petitive compared to RL-MRS which is designed to di- 
rectly handle Poisson noise. At low intensity regimes, our 
approach and RL-MRS are the two algorithms that give the 
best results. At high intensity, RL-TV performs very well, 
although RL-MRS, NaiveGauss and our approach are very 
close to it. NaiveGauss performs poorly at low intensity as 
it does not correspond to a degradation model with Pois- 
son noise. AnsGauss gives the worst results probably be- 
cause it does not handle the non-linearity of the degrada- 
tion model O after the VST. To assess the computational 
burden of the compared algorithms, Tab. [2 summarizes the 



execution times with an Intel PC Core 2 Duo 2GHz, 2Gb 
RAM. Except RL-MRS which is written in C++, all other 
algorithms were implemented in MATLAB. 

The same experiment as above was carried out with a 
microscopy image of the endothelial cell of the blood mi- 
crovessel walls; see Fig. [2 For the NaiveGauss, AnsGauss 
and our approach, the dictionary contained the wavelet 
orthogonal basis and the curvelet tight frame. The Ans- 
Gauss and the NaiveGauss results are spoiled by artifacts 
and suffer from a loss of photometry. RL-TV result shows 
a good restoration of small isolated details but with a dom- 
inating staircase-like artifacts. RL-MRS and our approach 
give very similar results although an extra-effort could be 
made to better restore tiny details. The quantitative mea- 
sures depicted in Fig. [3] confirm this qualitative discussion. 

Finally, we applied our algorithm on a real confocal 
microscopy image of neurons. Fig. Ua) depicts the ob- 
served imag^l using the GFP fluorescent protein. Fig.^b) 
shows the restored image using our algorithm with the or- 
thogonal wavelets. The images are shown in log- scale for 
visual purposes. We can notice that the background has 
been cleaned and some structures have reappeared. The 
spines are well restored and part of the dendritic tree is re- 
constructed, however some information can be lost (see 
tiny holes). This can be improved using more relevant 
transforms. 




Fig. 1. Deconvolution of a simulated neuron (Intensity ^ 30). (a) 
Original, (b) Blurred, (c) Blurred&noisy, (d) RL-TV, (e) NaiveGauss, (f) 
AnsGauss, (g) RL-MRS, (h) Our Algorithm. 

6. CONCLUSION 

In this paper, we presented a sparsity-based fast iterative 
thresholding deconvolution algorithm that take accounts of 
the presence of Poisson noise. Competitive results on con- 
focal microscopy images with state-of-the-art algorithms 
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Fig. 2. Deconvolution of the cell image (Intensity ^ 30). (a) Origi- 
nal, (b) Blurred, (c) Blurred&noisy, (d) RL-TV, (e) NaiveGauss, (f) Ans- 
Gauss, (g) RL-MRS, (h) Our Algorithm. 




Fig. 3. Mean t\ -error of all algorithms as a function of the intensity 
level for the deconvolution of the cell 

are shown. The combination of several transforms leads to 
some advantages, as we can easily adapt the dictionary to 
the kind of image to restore. The parameter A can be tricky 
to find, and we are developing a method helping to solve 
this issue. 
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Intensity regime 


Method 


<C 5 


<C 30 


<C 100 


<C 255 


Our method 


0.21 


0.76 


2.39 


4.79 


NaiveGauss 


0.41 


0.89 


2.31 


5.15 


AnsGauss 


2.07 


7.56 


21.25 


50.41 


RL-MRS 


0.21 


0.96 


2.2 


4.51 


RL-TV 


0.73 


1.47 


2.72 


4.28 



Table 1. Mean^i -error of all algorithms as a function of the intensity 
level for the deconvolution of the neuron phantom. 



Method 


Time (in s) 


Our method 


2.7 


NaiveGauss 


1.7 


AnsGauss 


1.7 


RL-MRS 


15 


RL-TV 


2.5 



Table 2. Execution time for the simulated neuron using the orthogonal 
wavelet transform 
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Fig. 4. Deconvolution of a real neuron, (a) Original noisy, (b) Restored 
with our algorithm 



